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We present atomistic molecular dynamics results for fully hydrated bilayers composed of ceramide 
NS-24:0, free fatty acid 24:0 and cholesterol, to address the effect of the different components 
in the stratum corneum (the outermost layer of skin) lipid matrix on its structural properties. 
Bilayers containing ceramide molecules show higher in-plane density and hence lower rate of passive 
transport compared to phospholipid bilayers. At physiological temperatures, for all composition 
ratios explored, the lipids are in a gel phase with ordered lipid tails. However, the large asymmetry 
in the lengths of the two tails of the ceramide molecule leads to a fluid like environment at the 
bilayer mid-plane. The lateral pressure profiles show large local variations across the bilayer for pure 
ceramide or any of the two component mixtures. Close to the skin composition ratio, the lateral 
pressure fluctuations are greatly suppressed, the ceramide tails from the two leaflets interdigitate 
significantly, the depression in local density at the inter- leaflet region is lowered, and the bilayer have 
lowered elastic moduli. This indicates that the observed composition ratio in the stratum corneum 
lipid layer is responsible for both the good barrier properties and the stability of the lipid structure 
against mechanical stresses. 



I. INTRODUCTION 

Stratum corneum (SC), the outer layer of the skin 
provides the main barrier against water loss and in- 
vasion by foreign pathogens. During their life-cycle, cells 
formed in the basal layer of the epidermis change their 
shape and composition. They progressively occupy outer 
layers of the epidermis, until they are peeled off from 
the outer layer. SC is often viewed as a bricks and mor- 
tar structure Q with corneocytes, the keratin filled non- 
viable disc like cells, arranged like bricks in a lipid mix- 
ture forming the mortar phase [i| . The three main com- 
ponents of the SC lipid matrix are a family of ceramide 
sphingolipids (CER), cholesterol (CHOL) and free fatty 
acid (FFA) 0, @ . Selective inhibition of any one of the 
ceramide Q, cholesterol Q or FFA Q is known to com- 
promise the barrier function of the skin. However, how 
the three components affect the lipid matrix properties 
at the molecular level is not known. To our knowledge, 
there have been only a few previous attempts at atom- 
istic modeling of the SC lipid layer. Most notably, Holtje 
et al. fioj simulated a two component mixture of fatty 
acid and cholesterol. Pandit and Scott [ll[ simulated a 
bilayer composed of symmetric CER NS 16:0 molecules. 
Notman et al. [i~2| used atomistic simulation to investi- 
gate 1 the effect of DMSO molecules on a hydrated bilayer 
composed of CER NS 24:0. But a systematic study of 
the effects of the three components is lacking and this 
study is specifically concerned with understanding the 
interplay of the three components which endows the skin 
with an almost contradictory combination of pliability 
and an extremely high penetration barrier. 

There are at least 9 different classes of ceramide found 
in human SC, with minor modifications in the head group 
region and the addition of an esterified fatty acid in the 
case of Ceramide 1. All the ceramides are conspicuous by 



having a large asymmetry in the length of the two tails 
and a large polydispersity in the fatty acid tail lengths 
[l3l ] . Similar polydispersity is found also in the length of 
the free fatty acids [lj]. Realistic representation of such 
a complex collection of molecules with atomic details is 
beyond current computational capabilities. Instead we 
choose just one representative ceramide, ceramide NS 
(also referred to as ceramide 2), with an asymmetric but 
monodisperse tail length. Ceramide NS is the most abun- 
dant species among the ceramide family. Its fatty acid 
tail is chosen to be 24:0, guided by the relative abundance 
of the different tail lengths in human SC [HI ; while its 
sphingosinc motif is chosen to have 18 carbons. Similarly 
we choose only FFA 24:0 because it is the most abundant 
free fatty acid found in SC lipid layer [l4j]. Fig Q] shows 
a skeletal representation of the molecules. 

Between the corneocytes, the lipid matrix shows regu- 
lar electron density variations, similar to lipid multilay- 
ers. This is not necessarily the only possible arrangement 
in the SC lipid matrix. In vitro exp eriments show the 
possibility of asymmetric leaflets [151 ] and multiple layer 
thicknesses, with indications that ceramide 1 connects 
different bilayers [l(J- The lipids are predominantly in 
a gel phase, possibly a single continuous gel phase [Tt} 
or with fluid regions [l^. In the skin, the corneocytes 
are arranged in clusters with the lipid matrix extend- 
ing through the full depth of the SC at intervals [l9l ]. 
These regions show a much lower permeation barrier than 
the layers between corneocytes [19( and it is not known 
whether or not the lipids there are arranged in a mul- 
tilayer structure. There is a hydration gradient across 
the SC, with the average water content around 30% by 
weight 0, [2l[ ■ How the water molecules are arranged in 
the SC is not completely known. 

With this uncertainty about the arrangement of the 
lipids in SC, and with the computational limitations im- 
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posed by detailed atomistic simulations, we study lipid 
mixtures arranged in symmetric bilayer structures in ex- 
cess water. An atomistic potential is chosen in this study 
to have a faithful representation of the lipids involved. In 
the absence of prior atomistic simulations, or, experimen- 
tal results on hydrated SC lipid bilayers, a coarse-grained 
description is not appropriate for this system. The water 
is required to stabilize the bilayer structure, and can be 
viewed as a replacement for the layering field imposed 
by the flat corneocytes. The lateral arrangement of the 
components and the properties of such a hydrated bi- 
layer will probably be representative of the lipid matrix 
confined between the corneocytes. 

The rest of the paper is organized as follows: In the 
section Simulation Details we outline the force field and 
protocols used for the simulations, and define the various 
properties measured. Next we provide the results for a 
wide range of compositions, concentrating mainly at tem- 
perature T=340K, which is much higher than the skin 
temperature. However, much of the results were found 
to remain qualitatively the same even at the skin tem- 
perature ~ 300K. We present the results for the higher 
temperature because we are more confident of the equi- 
libration at 340K than at skin temperature T~300K. Fi- 
nally we summarize the main findings and present an 
outlook for further studies on SC lipid structures. 



II. SIMULATION DETAILS 

A. Force field 

The interaction parameters used in the simulation are 
based on the united atom OPLS force field [HJ with mod- 
ifications for the nonpolar CH2/CH3 groups [23[ that 
accurately rep roduce experimental quantities for lipid 
molecules [24J. The polar hydrogens were included ex- 
plicitly and the partial charges for the headgroup atoms 
were selected to conform with molecules having similar 
structures and simulated using the same force field in 
the literature. The dihedral potentials in the hydrocar- 
bon tails were described by the Ryckaert-Bellemans term 

@. 

The force field used for the ceramide molecules has 
been used previously [l2|. The topologies of the polar 
part of the fatty acid and the cholesterol are the same 
as used injHj. The water is modeled with the SPC 
potential [2a ] . 



B. Simulations 

We used extended ensemble molecular dynamics sim- 
ulation at constant temperature and pressure (NPT) en- 
semble with GROMACS molecular dynamics package 
IH, ■ Temperature was controlled by Nose-Hoover 

thermostats [3l|, [32| coupled separately to the lipid and 
water molecules with time constant 5ps. Pressure was 
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FIG. 1: Skeletal representation of the atomic arrangements of 
the ceramide (CER NS 24:0), the free fatty acid (FFA 24:0) 
and cholesterol (CHOL). Partial charges used in the simula- 
tion are shown in brackets. Inset: (a) V-shaped and (b) hair- 
pin arrangements of the two tails of a CER molecule. [271 ] . 



controlled by an anisotropic Parrrinello-Rahman baro- 
stat (33l . [HI with time constant 5ps and compressibility 
4.5 x 10~ 5 /bar. The off-diagonal terms of the compress- 
ibility matrix were set to zero to preserve the orthogonal 
shape of the simulation box. Standard periodic boundary 
conditions were applied in all three directions. A simple 
group-based cut-off was used for calculating electrostatic 
interaction. The cut-off length was chosen to be 1.2 nm 
for both the Van der Waals and electrostatic potentials. 
Because the lipid molecules considered in this study do 
not fully ionize, the dominant electrostatic interaction is 
the dipole-dipole interaction, with the dipoles made up of 
the partial changes separated by fixed bond lengths. This 
is quite different from the case of phospholipids, where 
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fully ionized phosphate group introduces a large dipolc 
moment from well separated N~ and P + charges. As a 
specific example, the C25=026 atoms on CER molecule 
(see fig. 1) form a dipole from equal charges of magni- 
tude 0.28 e separated by a fixed distance of 0.123 nm. 
This creates a dipole moment of 1.65 D, which is roughly 
15 times smaller than the dipole moment at the head 
groups in DPPC bilayer [35| ■ Noting that the dipolar in- 
teraction falls of as 1/ r 3 and that for DPPC bilayer effect 
of electrostatic interaction with the periodic images be- 
come negligible for system sizes ~ 1024 lipid molecules 
[HI, we estimate that for the dipole moments involved 
in our system, the effect of electrostatic interaction with 
the periodic images become ncglibilc for systems contain- 
ing ~ 168 lipids. We have carried out simulations with 
pure CER molecules and 2:2:1 composition system both 
with group-based cut-off and particle mesh Ewald sum- 
mation (PME) for several system sizes. The results are 
presented in the appendix and validates the independence 
of the results presented on long range electrostatic inter- 
action and system size. Because the use of PME requires 
more than twice the computational cost as compared to 
using cut-off and it has no effect on the validity of the re- 
sults presented here, we restrict ourselves to group-based 
cut-off for handling electrostatic interaction in the simu- 
lations presented in the main paper. 



All bond lengths were constrained using the SHAKE 
algorithm [Hj]. The analytic SETTLE algorithm was 
used to handle the rigid SPC water molecules (37[- We 
used a timestep of 2 fs for T < 340 K and 1 fs for high 
temperature simulations. During the the production run, 
configurations at intervals of 0.5 ps were stored for fur- 
ther analysis. 

The SC lipid matrix (in mice) shows a pH of ~ 6 [38| . 
The ionization state of FFA is known to modify the struc- 
ture of CHOL-FFA (C16:0) mixtures [39j|. For ceramide- 
cholestcrol mixtures containing palmitic acid (C16:0) and 
oleic acid (C18:l), reported values of the pK a lie in the 
range of 6.2-7.3 [40|. Typically longer tailed FFA are ex- 
pected to have higher pK a ; e.g., for hexacosanoic acid 
(C26:0) in an egg-phosphatidylcholinc bilayer pK a ~ 7.4. 
This suggests that a fraction of the FFA should ionize at 
skin pH. However, the experimental evidence for the ef- 
fect of pH on SC lipid mixture is not clear. Surface force 
and AFM measurements on SC lipid mixtures containing 
free fatty acid seem to show no detectable effect in a wide 
range of pH (between 3.0 and 7.0) [4l|. The apparent in- 
sensitivity to pH for SC lipid mixture may be because the 
bilayers containing long tail ceramides remain in a dense 
gel phase and do not allow much structural freedom to 
the FFA molecules. With this background, we chose to 
simulate only un-ionized FFA - not having free charges 
allow fast calculation of the forces with cut-off as alluded 
before. 



C. Initial structures 

In all but one of the six different crystal structures 
formed under varying conditions, the sphingosine and 
the fatty acid tails of CER NS 24:0 arrange themselves 
with a large opening angle V-shaped structure p7| . Only 
low temperature crystallization from solution leads to a 
hairpin arrangement of the tails (inset of Fig. [!}. To 
ensure a hairpin arrangement of the ceramides, we start 
with a multilayer anhydride system with the ceramides 
placed in a slightly expanded hexagonal lattice with a 
hairpin structure at 300K, and relax the configuration 
with a 20 ns NPT MD run. A bilayer from the multi- 
layer was then placed in a larger box and the box was 
filled with water. We ran the system for 5 ns with the 
ccramide molecules fixed (to hydrate the bilayer prop- 
erly), and then for 2 ns with only the terminal methyl 
group on all the lipid tails frozen, so that the rest of the 
molecules reorient themselves to accommodate the water 
environment. This configuration was used as the initial 
configuration for the pure CER system. A further 10 ns 
NPT MD run was used at each temperature before the 
production runs. 
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TABLE I: Compositions (as molar ratios) and number of par- 
ticles simulated in this study. The pure FFA system had 9020 
water molecules. For all other compositions, there were 5250 
water molecules. 



For mixtures containing CER molecules, we transform 
the required number of CER molecules, chosen at random 
but maintaining the same composition in the two leaflets, 
to either CHOL or FFA. To convert to FFA molecules 
we simply separate the two tails of CER as two sepa- 
rate molecules and slowly grow the shorter chain derived 
from the sphingosine motif to the required length. FFA 
molecules are significantly more mobile as compared to 
CER molecules. Thus, although in the initial configu- 
ration the two FFA molecules generated from the same 
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CER molecule are next to each other, they do not intro- 
duce any artificial correlation after the equilibration step. 
To convert to CHOL we map certain atoms of the CER 
to CHOL. Repeated short NVT MD runs were performed 
while keeping the CHOL molecules frozen, and after each 
such short run, the atoms of the CHOL were displaced 
by a small amount until they reached the equilibrium po- 
sitions on the molecule. The mixtures were adequately 
equilibrated, typically by energy minimization followed 
by a series of short NVT simulations and finally with at 
least 10 ns of an NPT MD run. To find a typical re- 
laxation time, we quenched a CER bilayer from 360 K to 
300 K and found that the area/lipid and the bilayer thick- 
ness re-equilibrate to 300 K values in ~ 1 ns time-scale. 
For pure FFA or CHOL bilayers, we simulated multilay- 
ers and placed one bilayer in water in the same way as in 
the case of CER. Table Q] shows the number of molecules 
used, as well as the corresponding molar ratios, rounded 
to closest integer ratios. Some of the composition ratios 
explored in this study may not be achievable in exper- 
iments as hydrated bilayers. However, once prepared, 
the bilayers are kinetically stable and allow us to isolate 
the effects of the different components by studying the 
extreme compositions. In the rest of the paper the dif- 
ferent compositions arc referred to by the corresponding 
molar ratios. 



III. MEASURED QUANTITIES 

In this section we define the different quantities mea- 
sured from the simulations. For all composition ratios 
and temperatures, we use averages over 10 ns runs. We 
estimate the error bars for the measured quantities from 
the variance of intermediate averages over 2 ns windows. 

We define the bilayer normal direction to be the z di- 
rection. The orthogonal box shape along with the small 
system size ensures that this direction remains the same 
as the z direction of the simulation box. All z-dependent 
quantities reported are averaged over the lateral direction 
(x — y plane) in the entire simulation box. 



1. Structural properties 

At the molecular length scale, the water-lipid interface 
has a finite width. To assign an unique value to the bi- 
layer thickness 2d, we calculate the density of the water 
molecules as a function of z and define the position of the 
interface between the lipid and water as the z at which 
the water density decays to - of the bulk water density. 
We use this criteria, as opposed to the more usual cri- 
teria of the Gibbs dividing surface, because we consider 
multi-component lipids in this study and the different 
components have significantly different densities at the 
boundary. Concentrating on the approximately exponen- 
tially decaying water density gives a simple and unique 
method to assign a value to the bilayer thickness. The 



average lipid layer density p^ is computed by assuming 
that all the lipid mass is homogeneously distributed be- 
tween the two lipid- water interfaces. The local density of 
the lipid molecules depends on z with a minimum p™ m 
between the two leaflets. 

The asymmetry in the lengths of the two hydrocarbon 
tails of CER can lead to significant interdigitation. We 
define the following dimcnsionless overlap parameter as 
a quantitative measure of interdigitation, 



p t (z) x p b (z) 
[p t (z) + p b (z)Y 



(1) 



where pt(z) and Pb(z) are the densities at z for the CER 
from the top and bottom layers respectively. p v(z) = 1 
if half of the density at z is from the top layer and the 
other half is from the bottom layer CER molecules. If 
only the top or bottom layer of CER is at z, p v{z) = 0. 
Integrating over z, we define a single length scale X ov = 
Jo Pov{z)dz to compare the amount of interdigitation for 
different compositions. The integration is carried over 
the whole box, since if no CER is present, then p v{z) = 
and there is no contribution to X ov . 

The area compressibility ka of the bilayer is related to 
the area fluctuation in the NPT ensemble by [42j 



ka = k B T 



<A> 



< A 2 > - < A > 2 



(2) 



where kg is the Boltzmann's constant and T is the tem- 
perature. The angular brackets refer to averages over 
time for the area A and square of the area A 2 . 

The system sizes investigated in this work are too small 
to calculate the bending modulus by analyzing the un- 
dulation modes [H, |44[ ■ For fully saturated fluid mem- 
branes, the bending modulus n is approximately related 
to the area compressibility modulus ka and the thickness 
of the hydrocarbon tail region 2d c through 



K A {2d c f 



(3) 



where the constant c e is estimated to be 24 from a theory 
based on polymer brushes [45| . This simple prescription 
is found to work well for a number of different lipid bilay- 
ers Hll . In the present simulation, the head groups of 
the lipid molecules are much smaller than that of phos- 
pholipids and we use the bilayer thickness 2d as a good 
estimate of the hydrocarbon thickness 2d c . 



2. Tail order parameter 

The orientational (nematic) order of the tails is probed 
through an order parameter P 2 , defined by the largest 
eigenvalue of the second rank tensor 



Qa/3 = ( I -zUiaUifi 



-s af3 



(4) 
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where Ui a is the cartesian component a of some specific 
orientation vector on the lipid molecule i, and the < > 
denotes average over time and the lipid molecules. For 
fluid bilayers, Pi is related to the NMR deuterium or- 
der parameter Scd through the relation Pi = —2Scd 
[47I ] . The eigenvector corresponding to the largest eigen- 
value of Q defines the average orientation. Different car- 
bon atoms along the tail are expected to have different 
amount of ordering. To probe the local ordering, we de- 
fine a local orientation Pi(n) at carbon atom n on the 
lipid molecule i by the vector direction u ia (n) between 
carbon atoms n — 1 and n + 1. 



3. Local pressure 

The anisotropic arrangement of lipid molecules in a 
bilayer leads to an anisotropic local pressure profile, and 
hence local stress variations. While macroscopic pressure 
is a well defined quantity in molecular simulations, its 
microscopic description is not unique j4^|. We use the 
formalism of [4!| to define the local pressure tensor at a 
given height z as 



P(z) 



1 



E 



\i£slice i<j 

(5) 

Here, v$ is the velocity of particle i, is the force on 
particle i due to particle j and is the relative position 
vector of particle i from particle j . V s u ce is the volume 
of a thin slice with thickness Az and < > denotes an 
average over time. The first sum runs over the particles 
in the slice centered at z. The function / determines 
the contribution from the virial (the second sum) to the 
current slice. / is unity when both particles are in the 
current slice. If one or both of the particles are outside 
the slice (but the shortest distance between the two parti- 
cles goes through the slice) ; / is respectively chosen to be 
Az/\zi — Zj\ and Az/ (2|z, — Zj\). Because constraints in 
the simulation (SHAKE algorithm) transfer some of the 
kinetic contributions to the constraint potential, we need 
to consider both the sums in Eq. [5] explicitly. We cal- 
culate the local pressure profiles by starting with stored 
configurations separated by 2 ns and re-evolving the con- 
figurations for 200 ps. 

Of particular interest for lipid bilayers is the difference 
between the lateral and the normal pressure SP(z) = 
Plat(z) ~ Pzz(z), where P L at(z) = \ \P xx {z) + P yy (z)]. 
The surface tension 7 of the bilayer is related to SP(z) 
through [i| 
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5P(z)dz. 



(6) 



Anisotropic box fluctuations ensure that the bulk water 
is isotropic, i.e. the pressure tensor has equal diagonal 
components and zero off-diagonal components. This in 
turn leads to a zero average surface tension for the lipid 



bilayer. But locally 5P{z) goes through a number of 
maxima and minima. We define a microscopic stress by 



ep 



f 

2d 



dz 



(6P(z) - (SP(z)))' 



(7) 



Here, the integration is over the bilayer thickness and the 
angular brackets reflect time averages. 



IV. RESULTS AND DISCUSSION 



A. Pure Ceramide bilayers 





(a) 



(b) 




(c) 

FIG. 2: Snapshots of a CER2 bilayer at (a) 300K and (b) 
360K. While both the leaflets have similar tail ordering, pro- 
jection in two dimensions artificially accentuates order in one 
leaflet at the cost of showing less order in the other leaflet, 
(c) Time trace of two CER molecules. The perspective for the 
two molecules are chosen independently and each time frame 
has been shifted arbitrarily. 

Fig- HI shows snapshots of the CER bilayer at (a) 300K 
and at (b) 360K. The tails of the CER align along the 
z direction with strong ncmatic ordering. The terminal 
methyl groups at the bilayer midplane show a liquid-like 
disordered arrangement. With increasing temperature, 



the disordered region at the bilayer midplane gradually 
thickens. The tail order parameters at 300K and 360K 
are shown in Fig. [3l The terminal groups (smallest and 
largest atom indices) have small ordering. P2 also de- 
creases close to the head group (atom indices 16 and 25). 
The C-C bonds close to the head group orient at an angle 
to the vertical direction. Fig. 2] (a) shows a top view of 
one leaflet. The figure shows that the molecules arrange 
in a zigzag fashion with rows of molecules having the head 
group arrangement in alternate (orthogonal) directions. 
Thus the tail carbon bonds close to the head group point 
in one of these two orthogonal directions, reducing the 
overall average of P2. 
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FIG. 3: Tail order parameter of CER atoms as a function of 
atom index n (see Fig. [1]). 

In the well-ordered layer immediately beneath the head 
groups, the tails have fairly well-defined hexatic order 
with very few defects (Fig. [JJ)). Moreover, the orien- 
tations of the planes containing the lipid tails (Fig. 0] 
c ) suggests that the phase of this ordered layer is in- 
termediate between crystal and rotator phases (familiar 
from studies of alkanes of different lengths [I(|). While 
the two dimensional slices show strong sixfold order, the 
molecules undergo relatively fast slithering motion of the 
tails (Fig. [2]c). The coupling between the two leaflets 
is weak, allowing the two leaflets to oscillate about each 
other by more than the inter-lipid distance in roughly 
10 ns. Diffusion of the center of masses of the molecules 
in the x — y plane is much slower because they can only 
diffuse at the defect sites. 

Fig. [5] shows the density of the CER and water 
molecules across the bilayer at 300K (open symbols) and 
360K (closed symbols). The CER density shows a peak 
at water-lipid interface from the close packing of head 
group atoms. There is an almost constant density shoul- 
der from the dense packing of the hydrocarbon tail atoms, 
followed by a dip at the bilayer midplane, which corre- 
sponds to the amorphous inner layer due to the asym- 
metric ccramide tails. The arrow indicates where the 
water density falls below 1/e of the bulk water density at 




FIG. 4: (a) Top view of CER molecule head group arrange- 
ment on one leaflet at T = 340 K. The molecules are drawn 
in thick or thin lines respectively, depending on whether the 
tangent of the angle between the line joining the end-atoms 
of the head group and the x axis is positive or negative, (b) 
Cross section of CER tail groups a few atoms below the head 
groups, showing local hexatic order, (c) Side view of part 
of the lipid tails with (orthogonal distance regression) planes 
containing the tail atoms indicated by transparent boxes. The 
normals to these planes are indicated by solid rods on top of 
the boxes. In a crystalline configuration, the alignment of 
these planes is regular throughout the sample, while in the 
rotator phase, there are no correlations between the planes. 
In the present case, local correlation between the hydrocar- 
bon planes persist but allows slow rearrangement because the 
correlation is not perfect. 




FIG. 5: Density of CER atoms (circles), water molecules (tri- 
angles) for CER bilayer at 300K (filled symbols) and 360K 
(open symbol). The density of head group atoms at 300K is 
shown with a solid line. 
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300K, which is considered to be the lipid-water interface 
to calculate the bilayer thickness. 
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FIG. 6: Area/lipid (circles, left y-axis) and bilayer thickness 
(squares, right y-axis) of CER molecules as a function of tem- 
perature. 
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FIG. 7: Area compressibility (circles, left y-axis) and average 
density (squares, right y-axis) of CER bilayers as a function 
of temperature. 

Fig. [6] shows the area per lipid and bilayer thickness 
as a function of temperature. The area/lipid at 300K is 
0.389 nm 2 , which is close to the X-ray result ~ 0.4 nm 2 
for CER bilayer with molecules in hairpin arrangement 
[2t1 |. With increasing temperature, the area/lipid grad- 
ually increases, while the increasingly disordered termi- 
nal tail atoms make the thickness decrease. The average 
density and area compressibility (Fig. [7]) both decrease 
with temperature. Even at 360K, the area compressibil- 
ity is an order of magnitude larger (~ 4000 dyne/cm) 
than in most phospholipid fluid bilayers. The results 
do not show a sharp transition. Instead the disordered 



inter-leaflet region expands smoothly with temperature. 
The long chains involved in CER probably will lead to 
a broad gradual softening like a waxy material, instead 
of a sharp transition. Experimentally, main transition 
temperatures in ceramide systems have been found, for 
multilayer systems, to be of order 394 — 420 K (27|. One 
expects a lower temperature for a single hydrated mem- 
brane in solution. 



B. Mixed SC lipid bilayers 
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FIG. 8: Area compressibility (circles, left y-axis) and av- 
erage density (squares, right y-axis) of a mixed system, 
CER:CHOL:FFA=2:2:I, as a function of temperature. 

For all the composition ratios investigated the temper- 
ature behavior in the simulated range of temperatures is 
gradual and similar to that of CER system. (Fig.[H]shows 
results for the composition CER:CHOL:FFA=2:2:l). For 
presentation of the results in this section, we concentrate 
on 340K. A summary of the different measured quantities 
is included in Table HT1 

Fig. [9] shows snapshots (a) for pure CER, (b) equimo- 
lar CER-FFA, (c) equimolar CER-CHOL and (d) 2:2:1 
mixture of CER, CHOL and FFA. The tails of CER 
molecules (Fig. [5Ja,) retain substantial ordering. Long 
chain FFA molecules (Fig. [9}o) fall in registry with the 
CER molecules. The slightly longer length of the FFA 
molecules is accommodated by partially increasing the 
tail order in the CER molecules. The tails in the leaflets 
arrange in a slightly tilted orientation with respect to the 
layer normal. Thus, even though the amount of ordering 
is the same on both leaflets, the particular orientation of 
Fig. [5] highlights the ordering in the lower leaflets at the 
cost of obscuring the order in the upper leaflets. 

The head groups of the CHOL (Fig. [5};) mostly stay 
at the water-lipid boundary. The shorter length of the 
CHOL molecules squash the bilayer, with the central re- 
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composition 


2d 


PL 


pT d 




k a (10 a ) 


K (10 _ii ) 


ep 


(molar ratio) 


(nm) 


(g/ cc ) 


(g/cc) 


(nm) 


(dyn/cm) 


(erg) 


(bar) 


1:0:0 


5.66 


0.95 


0.66 


1.30 


4.9 


6.6 


660 


0:1:0 


3.16 


1.03 


0.91 


- 


4.7 


1.9 


870 


0:0:1 


6.33 


0.92 


0.51 




1.5 


2.6 


580 


7:1:0 


5.56 


0.95 


0.68 


0.93 


6.3 


8.1 


580 


3:1:0 


5.29 


0.95 


0.68 


1.16 


4.4 


5.2 


460 


2:1:0 


5.16 


0.95 


0.67 


1.42 


5.5 


6.1 


500 


1:1:0 


4.75 


0.96 


0.73 


1.87 


4.9 


4.6 


400 


7:0:1 


5.73 


0.95 


0.63 


0.90 


4.4 


6.0 


640 


3:0:1 


5.81 


0.95 


0.64 


1.10 


3.7 


5.2 


620 


2:0:1 


5.88 


0.95 


0.61 


0.90 


4.2 


6.0 


570 


1:0:1 


5.91 


0.94 


0.60 


1.42 


2.4 


3.5 


550 


1:1:1 


5.17 


0.94 


0.72 


2.20 


4.9 


5.5 


350 


1:2:1 


4.99 


0.96 


0.80 


2.52 


4.5 


4.7 


390 


2:1:1 


5.82 


0.93 


0.77 


3.33 


2.3 


3.3 


350 


5:5:1 


4.94 


0.95 


0.74 


1.48 


5.2 


5.3 


480 


2:2:1 


4.89 


0.94 


0.69 


2.50 


3.9 


3.9 


360 



TABLE II: Structural properties of the SC lipid bilayers for different ratios of CER:CHOL:FFA at 340K. 



gion comprising primarily CER tails. The ordering in the 
CER tails is decreased compared to pure CER bilayers 
and the tails from the two leaflets overlap strongly. 




(c) (d) 

FIG. 9: Snapshots of (a) pure CER, (b) equimolar CER-FFA 
mixture (1:0:1), (c) equimolar CER-CHOL mixture (1:1:0) 
and (d) 2:2:1 mixture of CER, CHOL and FFA at 340K. 
Only part of the water box is shown. CER, CHOL, and, FFA 
are colored as blue, orange, and green respectively. Water 
molecules are drawn with thin lines. 



The molecular arrangements in the ternary systems 
(Fig. (5J1) are somewhat intermediate to the binary 
CR2-FA and CR2-CHOL mixtures. Single chain FFA 
molecules are more flexible than CER molecules, which 
induces more FFA atoms in the midplane disordered 
phase and thus increases the tail order of CER molecules 
to some degree compared to CER-CHOL mixtures. 




2 4 -4 -2 a 

z (nm) z (nm) 



(c) (d) 

FIG. 10: Lipid densities along the bilayer normal directions 
at 340K for (a) CER (b) FFA (c) CHOL and (d) all lipid 
molecules considered together. Legends show molar ratio of 
CER:CHOL:FFA. 

To investigate the arrangement of the different com- 
ponents further, in Fig. [10] we plot the density of the 
lipid components along the z direction for different com- 
positions. Fig. [TUk shows the density profile of the CER 
atoms. Both pure CER and CER-FFA systems have a 
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n 



(nm) 



FIG. 11: Estimate of overlap \ v of the CER tails from the 
two leaflets (Eq.[TJ at 340K. 



high density near the head groups, followed by a con- 
stant density region due to the ordered tail atoms and 
then a region of lower density at the midplane, which 
covers almost 2nm. For CER-CHOL and three compo- 
nent mixtures, the density from the tail atoms of the 
CER molecules remains almost constant throughout the 
bilayer. In fact, for CER-CHOL mixtures, there is an in- 
crease in the local CER density at the bilayer midplane 
due to increased intcrdigitation of the molecules from the 
two leaflets. The plateau densities in Fig. [10] are fixed by 
the relative abundance of the different molecules in the 
mixtures. Fig. [TOb shows the density profile of the FFA 
atoms. For pure FFA and CER-FFA mixtures, the den- 
sity profiles are qualitatively similar to that of pure CER 
bilayers, except the region of low density intcr-leaflct 
space is narrower. In the 3-component mixtures, the den- 
sity minimum is replaced by a local density maximum be- 
cause of inter-leaflet overlaps. The density profile of the 
CHOL molecules (Fig. [TUb ) does not have the constant 
density hydrocarbon regions present for the FFA and 
CER molecules. The profile for the pure CHOL bilayer is 
not symmetric and the intcrlcaflet density quite high, sig- 
nifying significant transfer between the two leaflets dur- 
ing the simulation time scale. Fig. [TOH shows the total 
lipid densities. The effect of fatty acid is to increase 
the bilayer thickness and reduce the density by a small 
amount. CHOL, when present, reduces the bilayer thick- 
ness and increases the density at the tail region. 

The overlap of the lipids from the two leaflets (partial 
intcrdigitation) is expected to increase the intcr-leaflct 
friction and couple the dynamics of the two leaflets. 
Fig. [TT1 shows a surface plot of X ov , a measure of interdig- 
itation, as a function of composition. The three compo- 
nent mixture is characterized by large X ov . One of the SC 
lipids, not considered in this study, ceramide 1 (ceramidc 
EOS) contains an additional long-chain w-hydroxy acid 
(with number of carbon atoms > 30) linked to the fatty 
acid tail. This effectively makes the length of the long 



tail of ceramidc 1 nearly double to that of the other mem- 
bers of the ceramide family. The presence of ceramide 1 
will probably introduce significantly more interdigitation 
and inter-leaflet coupling. 
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FIG. 12: Local tail order parameter P2(n) as a function of 
atom index n on a CER molecule (Fig[T} for various compo- 
sitions at 340K. Errorbars are shown for the pure CER data 
(circles). 
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FIG. 13: Local tail order parameter Piiri) as a function of 
atom index n on a FFA molecule (FigfTJ for various compo- 
sitions at 340K. 

Fig. [T2] shows the tail order parameter of CER atoms 
at 340K for different composition ratios. The presence of 
CHOL reduces the nematic order, while FFA increases 
the order. The presence of either CER or CHOL re- 
duces the order of the FFA atoms compared to a pure 
FFA bilayer (Fig. [T3]) . This differs from phospholipid 
membranes, where the planar shape of cholesterol typ- 
ically increases the nematic order of the phospholipid 
tails. When mixed with the long ceramides, the shorter 



10 



cholesterol molecules tend to encourage a thinner mem- 
brane by disordering the longer ceramide tails so that 
they can fill the space around the cholesterol. This also 
accounts for the increased overlap of the CER tails upon 
adding cholesterol, as depicted in Fig. QT] 



The bending modulus k, calculated from the polymer 
brush theory [45| using ka and the bilayer thickness, be- 
haves in a similar fashion. Close to the skin composition 
CER:FFA:CHOL=2:2:l the bilayer becomes softer, with 
comparatively smaller ka and k. The absolute magni- 
tude of the clastic constants remain much higher than in 
the fluid phase of phospholipid bilayers. 
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FIG. 14: Area compressibility modulus ka of the bilayers at 
340K. 
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FIG. 15: Difference between the lateral and perpendicular 
components of the local pressure (8P = Plat — Pzz) as a 
function of distance from the bilayer mid-plane at 340K. (a) 
5P(z) for pure CER. (b) and (c) respectively show the effect 
of adding CHOL and FFA to CER bilayers. Because the 
local pressure profile is symmetric about the bilayer midplane, 
only one side of the data is reproduced for clarity, (d) SP 
for selected three component mixtures. For some data sets, 
smooth lines joining the points are drawn as a guide to eye. 
Legends show the ratio of CER:CHOL:FFA. 

Fig. [TJ] shows a surface plot of the area compressibil- 
ity ka on a ternary diagram of the three components. 
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FIG. 16: Local stress Ip (Eq.[7]l in the bilayers associated with 
the difference in the lateral and the perpendicular components 
of the pressure tensor. 2:2:1 composition is indicated with a 
vertical line. 

Fig. rT5] shows the difference between the lateral and 
the normal pressure 5P(z) along the bilayer normal di- 
rection. Inside the bilayers, pure CER (Fig. \15k ) shows 
large variations in SP(z) as a function of z. Similar but 
less pronounced variations in SP{z) are observed in phos- 
pholipid membranes. For example, Lindahl and Olle (44| 
found a maximum variation of |(5P(z)| < 500 bar for a 
DPPC bilayer at 323K, which is 5 times smaller than 
the present case. The addition of CHOL (Fig. [TSb) in- 
troduces additional peaks in SP(z). FFA reduces the 
peak height marginally (Fig. [15b). The three component 
mixtures (|15H ) show much less pronounced variations in 
SP(z) than the pure ceramide or any of the 2 component 
mixtures. A large SP can be interpreted as a local mo- 
ment acting on the molecules, making a deviation from 
the flat interface of bilayer more likely. The magnitude of 
SP also can be viewed as a local stress ip. Fig. [TBI shows 
a ternary plot of ip defined through Eq. [7] The vertical 
line shows the position of composition 2:2:1, where ip 
shows a minimum. To gain some insight into the energy 
scale, we note that 10 3 bar ~ 0.7/csT per methyl group, 
while the trans-gauche energy difference is ~ 2ksT per 
dihedral bond. Thus, in the absence of prohibitive en- 
ergy barriers, the system can introduce gauche defects to 
reduce the local variation of SP. 
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V. CONCLUSIONS 



We have presented molecular dynamics simulation re- 
sults for various composition ratios of the three main 
constituent components of stratum corneum lipid layers, 
namely, CER NS 24:0, FFA 24:0 and CHOL. The long 
asymmetric tails of the CER molecules form a dense bi- 
layer phase in water, with considerable interdigitation of 
the tails from the two leaflets and strong nematic order. 
CHOL, being rigid and smaller in length compared to 
CER, acts as a molecular clamp, squashing the bilayer 
and increasing the in-plane density. The resulting struc- 
ture shows large lateral stress variations, which is relieved 
by the presence of FFA. The three component mixtures 
are characterized by a higher density and comparatively 
smaller area compressibilities and bending moduli. 

The lipids in real skin stratum corneum have large 
polydispersity in tail length and have different head 
groups. This work does not address the effect of poly- 
dispersity. We also study bilayers in excess water, which 
is quite different from limited water environment in the 
skin. The strong local pressure fluctuations in the bi- 
layer suggests that the bilayer may not be the most stable 
structure and it will be interesting to simulate the lipids 
in limited water conditions. The hydrocarbon densities 
in the bilayers were found to be quite large compared to 
phospholipids. This is probably responsible for three or- 
ders of magnitude smaller permeability of water through 
stratum corneum as compared to plasma membranes. 



APPENDIX A: EQUILIBRATION TIME 

To probe the typical time it takes for the system to 
adapt the equilibrium conformation for a given temper- 
ature, we used an equilibrated configuration of 128 CER 
lipids (1:0:0 composition) at 360 K and monitored the 
instantaneous arca/lipid as a function of time (small 
black circles in Fig [17)) after abruptly changing the de- 
sired thermostat temperature to 300 K. We also show the 
area/lipid calculated from separate runs at 360 K and at 
300 K as shaded big circles with error bars at the two 
extreme ends of the data. As can be seen from the plot, 
the area/lipid adapts to the new temperature with a time 
scale of around 1 ns (position of the arrow) . In the paper 
we use equilibration times of 10 ns for each 10K temper- 
ature shift - which gives ample time for the molecules to 
reorganize themselves into the equilibrium conformation 
corresponding to the set temperature. 



APPENDIX B: EFFECT OF LONG RANGE 
ELECTROSTATICS AND SYSTEM SIZE 

To explore the dependence of measured quantities on 
system size and run times, we have carried out some 
additional simulations on pure CER bilayers (at 300 K) 
and bilayers with 2:2:1 composition ratio of CER, CHOL, 
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FIG. 17: Equilibration of area/lipid for a CER bilayer after 
a sudden change in thermostat temperature from 360 K to 
300 K. The big circles with error-bars are the values obtained 
from different runs at 360K and 300K. 



FFA (at 340 K). These simulations use much larger num- 
ber of lipid molecules and much longer run times than the 
simulations reported in the paper. For these simulations, 
we use both the group-based cut-off and PME methods 
to probe the effect of long-range electrostatics. In simula- 
tion with PME, the grid spacing was chosen to be 0.1 nm 
and fourth order polynomial interpolation was used. A 
summary of the number of molecules used and timescales 
probed is given in Tabic. [TTT1 To get the initial configura- 
tion for a larger system, we use the already equilibrated 
configuration at the required temperature from the simu- 
lations presented in the main paper and replicate it in the 
x and y directions. To reduce disk space requirements, in 
most cases of these set of simulations, we save configura- 
tions only at intervals of 0.2 ns (as opposed to 0.5 ps in 
the main paper) . Thus we expect larger statistical errors 
in the data reported in this appendix as compared to the 
main paper. 



Composition 


Number of molecules 


Run time (ns) 


CER:CHOL:FFA 














(molar ratio) 


CER 


CHOL 


FFA 


SOL 


Cutoff 


PME 


2:2:1 


56 


56 


32 


5250 


100 


100 


(340 K) 


224 


224 


128 


21000 


70 


40 




504 


504 


288 


47250 


50 


16 


1:0:0 


128 






5250 


30 


30 


(300 K) 


512 






21000 


30 


16 




1152 






47250 


20 


10 



TABLE III: Number of molecules and simulation times used 
to probe effect of long range electrostatics and system size. 
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In Tabic HVl wc report the energies (normalized by the 
number of lipid molecules in the system) for 2:2:1 system 
calculated separately with group based cut-off and parti- 
cle mesh Ewald summation (PME). The two important 
observations from the energy values are that (i) the differ- 
ence between the group based cut-off and PME schemes 
in the total energy is about 2% and (ii) the normalized 
energies arc independent of the system sizes. Both of 
these two observations show that the interactions with 
the periodic images do not contribute significantly in the 
electrostatic energy. Table IIVI also shows the area com- 
pressibility and the average bilayer thickness calculated 
from these simulations. The effect of the system size and 
handling of electrostatics do not affect the values signifi- 
cantly. 



Number of 


Coulomb 


Coulomb 


Total 


KA T 


2 d 


lipids 


(SR)* 


(LR)* 


energy* 




(nm) 


144 


cut-off 


-1.7675 




-1.1310 


2.6 (2) 


4.90 (2) 




PME 


-1.6363 


-0.1101 


-1.1114 


2.0 (3) 


4.95 (2) 


576 


cut-off 


-1.7676 




-1.1308 


2.5 (5) 


4.97 (2) 




PME 


-1.6364 


-0.1101 


-1.1105 


3.6 (3) 


4.91 (2) 


1296 


cut-off 


-1.7676 




-1.1302 


2.3 (5) 


4.93 (2) 




PME 


-1.6364 


-0.1102 


-1.1104 


2.9 (4) 


4.90 (1) 



TABLE IV: Short range (SR) and long range (LR) electro- 
static energies, total energy, area compressibility (ka) and 
bilayer thickness (2d) for bilayers containing 2:2:1 molar ratio 
of CER2, CHOL and FFA at 340K. Statistical errors of the 
mean values are indicated in the brackets as the uncertainty 
on the last digit. (* Energies are in MJ/mol and are normal- 
ized by the number of lipid molecules in the system, t ka is 
in 10 3 dyn/cm.) 

Table [V] shows the effect of system size and method 
of calculating electrostatic interaction for pure CER bi- 
layers at 300K. The conclusions about insensitivity of 
the results on both the system size and the long range 
electrostatic interaction in the case of 2:2:1 lipid mixture 



remains equally valid for the case of CER bilayers. 



Number of 


Coulomb 


Coulomb 


Total 


K ,t 

i" A 


2 d 


lipids 


(SR)* 


(LR)* 


energy* 




(nm) 


128 


cut-off 


-2.2045 




-1.5516 


7.0 (6) 


5.67 (1) 




PME 


-1.9974 


-0.1755 


-1.5272 


6.1 (6) 


5.68 (1) 


512 


cut-off 


-2.1996 




-1.5520 


6.2 (9) 


5.70 (1) 




PME 


-1.9973 


-0.1754 


-1.5269 


7.4 (9) 


5.70 (1) 


1152 


cut-off 


-2.1994 




-1.5516 


5.8 (7) 


5.72 (1) 




PME 


-1.9946 


-0.1753 


-1.5245 


7.4 (9) 


5.72 (1) 



TABLE V: Short range (SR) and long range (LR) electro- 
static energies, total energy, area compressibility (ka) and 
bilayer thickness (2c?) for pure CER bilayers at 300K. Statis- 
tical errors of the mean values are indicated in the brackets as 
the uncertainty on the last digit. (* Energies are in MJ/mol 
and are normalized by the number of lipid molecules in the 
system. ' ka is in 10 3 dyn/cm.) 

For the smallest system sizes used in this study, typ- 
ical CPU requirements for PME based calculations are 
more than twice as compared to the CPU requirement 
for group-based cut-off. Hence, for the simulation re- 
sults presented in the main paper, we confine ourselves 
to group-based cut-off scheme only. 
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